Section 9 Model Parameterization

These sections are not detailed enough to warrant their own dedicated chapter.

9.1 Parameter Changes

These sections mostly involve changing a few parameter values and are therefore in small sub sections without R-code. Some of this might change later.

9.1.1 Channel parameters revised

This an ongoing unresolved issue #49.

So far nothing has been changed.

9.1.2 Crop parameters verified

The crop database has been updated and stored in the following file:

"model_data/cs10_setup/swat_input_mod/plants.plt"

Discussions to this topic can be found in issue #27.

From the whole database, we only use crops defined in our management files:

Crop management: oats, barl, wwht, swht, pota, fesc

Generic: frst, fesc and others which are minor

The parameters for these crops have been updated to reflect the colder growing conditions.

9.1.3 Soil physical parameters in final form

This has been mentioned in issue #28, and has been closed without discussion. Seems like the parameters are in final form, but no further documentation exists at this point.

9.1.4 Soil chemical parameters in final form

This step will be verified in issue #46. It has been discussed in #19.

The following changes have been implemented, but are subject to change:

Changes made to soilnut1 of nutrients.sol
Parameter Default CS10
lab_p 5 20
nitrate 7 8.5
inorgp 3.5 35.1
watersol_p 0.15 0.4
nutrients_sol <- readLines("model_data/cs10_setup/swat_input/nutrients.sol")

header <- nutrients_sol[1]

nutrients_sol_df <- read.table("model_data/cs10_setup/swat_input/nutrients.sol", header = T, skip = 1, sep = "", fill = T, as.is = T, colClasses = "character")

nutrients_sol_df$lab_p <- 20
nutrients_sol_df$nitrate <- 8.5
nutrients_sol_df$inorgp <- 35.1
nutrients_sol_df$watersol_p <- 0.4

write.table(
  nutrients_sol_df,
  "model_data/cs10_setup/swat_input_mod/nutrients.sol",
  quote = F,
  sep = "   ",
  row.names = F, 
)

header_new <- paste(header, "and modified by the CS10 workflow on", Sys.time())

nutrients_sol <- readLines("model_data/cs10_setup/swat_input_mod/nutrients.sol")

writeLines(c(header_new, nutrients_sol), "model_data/cs10_setup/swat_input_mod/nutrients.sol")

Testing if the setup still works:

simout <- test_swat_mod()
cat(tail(simout), sep = "\n")

9.1.5 Impoundment parameters defined

Has been postponed by #20, will be dealt with in #54

9.1.6 Water diversions defined

Deemed as not relevant for this catchment

9.1.7 Point sources parameters added

Is an ongoing issue in #21

9.1.8 Tile drainage parameters defined

Is an ongoing issue in #53.

old_path <- "model_data/cs10_setup/swat_input/tiledrain.str"
new_path <-  "model_data/cs10_setup/swat_input_mod/tiledrain.str"

tiledrain_str <- readLines(old_path)

header <- tiledrain_str[1]

header <- paste(header, "then modified by CS10 workflow at", Sys.time())

tiledrain_df <-
  read.table(
    old_path,
    sep = "",
    header = T,
    skip = 1,
    fill = T,
    colClasses = "character"
  )

tiledrain_df$dp <- 800 # from 1000
tiledrain_df$t_fc <- 12 # from 24
tiledrain_df$lag <- 30 # from 96
tiledrain_df$rad <- 200 # from 100
tiledrain_df$dist <- 8000 # from 30
tiledrain_df$drain <- 40 # from 10
tiledrain_df$pump <- 0 # from 1
tiledrain_df$lat_ksat <- 2 # from 2 (no change)

write.table(header, file = new_path, quote = F, col.names = F, row.names = F)
write.table(
  x = tiledrain_df,
  file = new_path,
  sep = "   ",
  quote = F,
  append = T,
  row.names = F
)
## Warning in write.table(x = tiledrain_df, file = new_path, sep = " ", quote = F,
## : appending column names to file

Testing if the setup still works:

simout <- test_swat_mod()
cat(tail(simout), sep = "\n")

9.1.9 Atmospheric deposition defined

Has been done in section REF

9.1.10 Additional settings verified

Refers to chapter 3.11 in The Protocol and files parameters.bsn codes.bsn Has been discussed in #22

The following changes have been made:

Changes made to codes.bsn
Parameter Default CS10
pet 1 2
rte_cha 0 1
cn 0 1
tiledrain 0 1
soil_p 0 1
atmo_dep a m
old_path <- "model_data/cs10_setup/swat_input/codes.bsn"
new_path <- "model_data/cs10_setup/swat_input_mod/codes.bsn"

codes_bsn <- readLines(old_path)

header <- codes_bsn[1]

header <- paste(header, "then modified by CS10 workflow on", Sys.time())

codes_df <- read.table(old_path, skip = 1, header = T, sep = "", 
                       colClasses = "character")

PET method (pet) The OPTAIN project recommends the PET calculation of Pennman-Monteith, which is code 2.

codes_df$pet <- 2 

Channel routing network (rte_cha) Recommended is to start with Muskingum (code 1) and apply variable storage method if it causes problems.

codes_df$rte_cha <- 1

Stream water quality (wq_cha) Recommended to test both for OPTAIN (1/0)

Daily curve number calculation (cn) Code 1 is recommended. (Plant ET)

codes_df$cn <- 1

OPTAIN recommends new version

codes_df$soil_p <- 1

Lapse rate is being tested in REF

Plant growth stress is being testing in REF

Tile drainage

This should be set to 1, but that causes a crash. We will fix this bug in a dedicated secton (REF)

codes_df$tiledrain <- 0

Atmosdep was done in REF and is annual

codes_df$atmo_dep <- "y"

Writing the changes

write.table(header, file = new_path, col.names = F, row.names = F, quote = F) 


write.table(
  codes_df,
  file = new_path,
  col.names = T,
  append = T,
  quote = F,
  sep = "   ",
  row.names = F
)
## Warning in write.table(codes_df, file = new_path, col.names = T, append = T, :
## appending column names to file

Testing if the setup still works:

simout <- test_swat_mod()
cat(tail(simout), sep = "\n")

9.2 Fixing Tiledrain

run_this_chapter = FALSE

Problem: tiledrain=1 in codes.bsn causes crash.

Reason: Some HRUs have drains, but are on a shallow soil type. This causes the drains to be below the soil, leading to a (nondescript) error.

Likely source of error: Misalignment of the soil and land use map.

Discussion of this issue can be found in #53

require(mapview)
require(sf)
require(stringr)
require(dplyr)
require(purrr)
require(DT)

9.2.1 Finding the problem HRUs

Our drains are set to 80cm depth. Which of our soils are less than that?

usersoil <- readr::read_csv("model_data/input/soil/UserSoil_Krakstad.csv", show_col_types = F)

shallow_soils <- usersoil$SNAM[which(usersoil$SOL_ZMX<800)]

print(shallow_soils)
##  [1] "Mountain"     "CMlen-dy"     "Antropogenic" "Sea_dep_thin" "Humic"       
##  [6] "LVlen-sl"     "End_moraine"  "RGlen-dy-ar"  "RGlep-dy"     "STlen-dy"    
## [11] "STlen-lv-sl"  "STlen-rt-sl"  "STlen-um-sl"

Which of our HRUs use these soils, and are drained?

hru_data <- read.table( "model_data/cs10_setup/swat_input/hru-data.hru",
                        skip = 1, header = T, sep = "")

shallow_hrus <- which(hru_data$soil %in% shallow_soils)

drained_hrus <- grepl(x = hru_data$lu_mgt, pattern = "_drn") %>% which()

# Intersection of HRUs which are both shallow and drained
problem_hrus <- base::intersect(shallow_hrus, drained_hrus)

length(problem_hrus)
## [1] 16

Just 16 HRUs have an issue. Lets get some info on these 16:

problem_hru_data <- hru_data[problem_hrus,]

datatable(problem_hru_data)

We need to link the soil type to the hru vector file

hru_shp <- read_sf("model_data/cs10_setup/optain-cs10/data/vector/hru.shp")
hru_soil <- hru_data %>% select(soil, name)
hru_shp <- left_join(hru_shp, hru_soil, by = "name")

And separate out our problematic HRUs and affected LUs

problem_shps <- hru_shp %>% filter(name %in% problem_hru_data$name)
problem_hru_ids <- problem_hru_data$lu_mgt %>% str_remove("_drn_lum")
problem_lus <- hru_shp %>% filter(type %in% problem_hru_ids)

Next we remove the problem HRUs and affect LUs from our hru_shp (this is for the upocoming map, to remove overlaps).

hru_shp <- hru_shp %>% filter((hru_shp$name %in% problem_lus$name) == FALSE)
hru_shp <- hru_shp %>% filter((hru_shp$name %in% problem_shps$name) == FALSE)

We also remove the problem HRUs from the affected LUs

problem_lus <- problem_lus %>% filter((problem_lus$name %in% problem_shps$name) == FALSE)

Data processing done, lets have a look:

hru_map <-mapview(
  hru_shp,
  map.types = "Esri.WorldImagery",
  color = "green",
  lwd = 1,
  legend = FALSE,
  zcol = "soil",
  alpha.region = 0
)

problem_map <-
  mapview(
    problem_shps,
    color = "red",
    lwd = 3,
    legend = FALSE,
    zcol = "name",
    layer.name = "Problem HRUs",
    alpha.region = 0
  )

problem_lus_map <- 
  mapview(
    problem_lus,
    color = "orange",
    lwd = 1,
    legend = FALSE,
    zcol = "name",
    layer.name = "Problem LUs",
    alpha.region = 0
  )

full_map <- hru_map + problem_lus_map+ problem_map

full_map

Resolution to this problem has been discussed in issue #53

We will extract the hru and landuse file to modify further.

lum_fixed <- read.table("model_data/cs10_setup/swat_input_mod/landuse.lum", skip = 1, sep = "", header = T) %>% tibble()

hru_fixed <- read.table("model_data/cs10_setup/swat_input/hru-data.hru", skip = 1, sep = "", header = T) %>% tibble()

9.2.2 Fixing the problem HRUs

hru0218 with landuse a_008f_2_drn_lum – Remove the drains since it is an isolated LU

lum_fixed$tile[which(lum_fixed$name=="a_008f_2_drn_lum")] = "null"

hru3637 with landuse a_148f_1_drn_lum – Remove the drains since it is an isolated LU.

lum_fixed$tile[which(lum_fixed$name=="a_148f_1_drn_lum")] = "null"

hru3686 of landuse a_115f_1_drn_lum has the wrong soil type, changing it to "STlv-sl” of its connected field.

hru_fixed$soil[which(hru_fixed$name=="hru3686")] = "STlv-sl"

hru4880 of landuse a_182f_5_drn_lum has the wrong soil type for the land use. Changing to the “ARdy” of neighboring fields

hru_fixed$soil[which(hru_fixed$name=="hru4880")] = "ARdy"

hru5033 with isolated landuse a_112f_1_drn_lum can have its drains removed.

lum_fixed$tile[which(lum_fixed$name=="a_112f_1_drn_lum")] = "null"

hru0692 of isolated land use a_037f_1_drn_lum can safely has its drains removed

lum_fixed$tile[which(lum_fixed$name=="a_037f_1_drn_lum")] = "null"

hru0704 has the wrong soil type for its landuse a_037f_2_drn_lum, which should be STlv-sl instead of Sea_dep_thin

hru_fixed$soil[which(hru_fixed$name=="hru0704")] = "STlv-sl"

hru3314 of landuse a_060f_2_drn_lum has the wrong soil type. It should be STlv-sl, not STlen-rt-sl.

The same is true for hru3376.

hru_fixed$soil[which(hru_fixed$name=="hru3314")] = "STlv-sl"
hru_fixed$soil[which(hru_fixed$name=="hru3376")] = "STlv-sl"

hru2784 (and the within hru3367) is a sizable field with landuse a_072f_1_drn_lum and soil RGlen-dy-ar. Tough call here, but best bet is probably to assign soil of the other major field it shares a landuse with (STlv-sl)

hru_fixed$soil[which(hru_fixed$name=="hru2784")] = "STlv-sl"
hru_fixed$soil[which(hru_fixed$name=="hru3367")] = "STlv-sl"

hru3465 of LU a_146f_2_drn_lum has the wrong soil for the greater field. Changing it from RGlep-dy to STrt-sl

hru_fixed$soil[which(hru_fixed$name=="hru3465")] = "STrt-sl"

hru5214 of isolated land use a_020f_4_drn_lum can safely have its drains removed.

lum_fixed$tile[which(lum_fixed$name=="a_020f_4_drn_lum")] = "null"

hru5861 of landuse a_210_f2_drn_lum does not share the same soil (end moraine) as the rest of its field.

hru_fixed$soil[which(hru_fixed$name=="hru5861")] = "STlv-sl"

hru5926 of landuse a_211f_1_drn_lum has the wrong soil type (Sea-dep-thin), changing to that of the greater field (STrt-sl)

hru_fixed$soil[which(hru_fixed$name=="hru5926")] = "STlv-sl"

hru2633 of isolated LU a_131f_1_drn_lum can safely have its drains removed.

lum_fixed$tile[which(lum_fixed$name=="a_131f_1_drn_lum")] = "null"

9.2.3 Writing modified changes

Removing the appended columns, and duplicates from the left-join, adding header and writing to modified input files directory.

lu_header <- paste("lu table after fixing up by OPTAIN-CS10 workflow on", Sys.Date())

write.table(lu_header, "model_data/cs10_setup/swat_input_mod/landuse.lum", quote = F, row.names = F, col.names = F )

write.table(lum_fixed, file = "model_data/cs10_setup/swat_input_mod/landuse.lum", sep = "\t", quote = F, append = T, row.names = F, col.names = T)
## Warning in write.table(lum_fixed, file =
## "model_data/cs10_setup/swat_input_mod/landuse.lum", : appending column names to
## file

Six drains have been removed.

hru_header <- paste("header header header, delete me and regret it forever! fixed soil types by cs10 worflow on", Sys.Date())

write.table(hru_header, "model_data/cs10_setup/swat_input_mod/hru-data.hru", quote = F, row.names = F, col.names = F )

write.table(hru_fixed, file = "model_data/cs10_setup/swat_input_mod/hru-data.hru", sep = "\t", quote = F, append = T, row.names = F, col.names = T)
## Warning in write.table(hru_fixed, file =
## "model_data/cs10_setup/swat_input_mod/hru-data.hru", : appending column names
## to file

Ten soil types have been changed.

9.2.4 Re-enabling tiledrain

Load from modified input files

old_path <- "model_data/cs10_setup/swat_input_mod/codes.bsn"
new_path <- "model_data/cs10_setup/swat_input_mod/codes.bsn"

codes_bsn <- readLines(old_path)

header <- codes_bsn[1]

header <- paste("modified by CS10 workflow on", Sys.Date())

codes_df <- read.table(old_path, skip = 1, header = T, sep = "", 
                       colClasses = "character")

Re-enable by setting to 1

codes_df$tiledrain <- 1

Writing to modified input files directory.

write.table(header, file = new_path, col.names = F, row.names = F, quote = F) 

write.table(
  codes_df,
  file = new_path,
  col.names = T,
  append = T,
  quote = F,
  sep = "   ",
  row.names = F
)
## Warning in write.table(codes_df, file = new_path, col.names = T, append = T, :
## appending column names to file

9.2.5 Testing SWAT+

Now we can see if the model runs

source("model_data/code/test_swat_mod.R")

simout <- test_swat_mod()
cat(tail(simout), sep = "\n")

Success.